arXiv:1509.03162v4 [q-bio.NC] 2 Mar 2017 


Fundamental activity constraints lead to 
specific interpretations of the connectome 

Jannis Schuecker^®, Maximilian Schmidt Sacha J. van Albada^, Markus 

Diesmann^’^’^, and Moritz Helias^’^ 


^Institute of Neuroscience and Medicine (INM-6) and Institute for Advanced Simulation (IAS-6) 
and JARA BRAIN Institute I, Jiilich Research Centre, 52428 Jiilich, Germany 

^Department of Psychiatry, Psychotherapy and Psychosomatics, Medical Faculty, RWTH Aachen 
University, 52062 Aachen, Germany 

^Department of Physics, Faculty 1, RWTH Aachen University, 52062 Aachen, Germany 

®Both authors contributed equally to this work. 


Correspondence to: Jannis Schuecker, Wilhelm-Johnen-Strafie, 52425 Jiilich 
j .schuecker@fz-juelich.de 


Activity constraints on the connectome 


Schuecker, Schmidt et al. 


Abstract 

The continuous integration of experimental data into coherent models of the brain is an increasing chal¬ 
lenge of modern neuroscience. Such models provide a bridge between structure and activity, and identify 
the mechanisms giving rise to experimental observations. Nevertheless, structurally realistic network 
models of spiking neurons are necessarily underconstrained even if experimental data on brain connec¬ 
tivity are incorporated to the best of our knowledge. Guided by physiological observations, any model 
must therefore explore the parameter ranges within the uncertainty of the data. Based on simulation 
results alone, however, the mechanisms underlying stable and physiologically realistic activity often re¬ 
main obscure. We here employ a mean-field reduction of the dynamics, which allows us to include 
activity constraints into the process of model construction. We shape the phase space of a multi-scale 
network model of the vision-related areas of macaque cortex by systematically refining its connectivity. 
Fundamental constraints on the activity, i.e., prohibiting quiescence and requiring global stability, prove 
sufficient to obtain realistic layer- and area-specific activity. Only small adaptations of the structure are 
required, showing that the network operates close to an instability. The procedure identifies components 
of the network critical to its collective dynamics and creates hypotheses for structural data and future 
experiments. The method can be applied to networks involving any neuron model with a known gain 
function. 


Author summary 

The connectome describes the wiring patterns of the neurons in the brain, which form the substrate 
guiding activity through the network. The influence of its constituents on the dynamics is a central topic 
in systems neuroscience. We here investigate the critical role of specific structural links between neuronal 
populations for the global stability of cortex and elucidate the relation between anatomical structure and 
experimentally observed activity. Our novel framework enables the evaluation of the rapidly growing 
body of connectivity data on the basis of fundamental constraints on brain activity and the combination 
of anatomical and physiological data to form a consistent picture of cortical networks. 

Introduction 

The neural wiring diagram, the connectome, is gradually being filled through classical techniques com¬ 
bined with innovative quantitative analyses (Markov et al., 2014a,b) and new technologies (Wedeen et ah, 
2008; Axer et al., 2011). The connectivity between neurons is considered to shape resting-state and task- 
related collective activity (Cole et ah, 2014; Shen et ah, 2015). For simple networks, clear relationships 
with activity are known analytically, e.g., a dynamic balance between excitatory and inhibitory inputs 
in inhibition-dominated random networks leads to an asynchronous and irregular state (van Vreeswijk & 
Sompolinsky, 1996; Amit & Brunei, 1997a; Tetzlaff et ah, 2012). Structures and mechanisms underlying 
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Figure 1 The integrative loop between experiment, model, and theory. Black arrows represent the 
classical forward modeling approach: Experimental anatomical data are integrated into a model structure, 
which gives rise to the activity via simulation. The model activity is compared with experimental functional 
data. The usual case of disagreement leads to the need to change the model definition. By experience and 
intuition the researcher identifies the parameters to be changed, proceeding in an iterative fashion. Once the 
model agrees well with a number of experimental findings, it may also generate predictions on the outcome of 
new experiments. Red arrows symbolize our new approach: informed by functional data, an analytical theory 
of the model identifies critical connections to be modified, refining the integration of data into a consistent 
model structure and generating predictions for anatomical studies. 


large-scale interactions have been identified by means of dynamical models (Deco & Jirsa, 2012; Cabral 
et ah, 2014) and graph theory (Markov et ah, 2013; Goulas et ah, 2014). Furthermore, the impact of 
local network structure on spike-time correlations is known in some detail (Ostojic et ah, 2009; Pernice 
et ah, 2011; Trousdale et ah, 2012). Under certain conditions, there is a one-to-one mapping between cor¬ 
relations in neuronal network activity and effective connectivity, a measure that depends on the network 
structure and on its activity (van Albada et ah, 2015). In conclusion, anatomy does not uniquely deter¬ 
mine dynamics, but dynamical observations help constrain the underlying anatomy. Despite advances in 
understanding simple networks, a complete picture of the relationship between structure and dynamics 
remains elusive. 

Fig. 1 visualizes the integrative loop between experiment, model, and theory that guides the inves¬ 
tigation of structure and dynamics. In the traditional forward modeling approach, structural data from 
experimental studies determine the connectivity between model neurons. Combined with the specification 
of the single-neuron dynamics and synaptic interactions, simulations yield a certain network dynamics. 
There is a fundamental problem with this approach. 

Despite the impressive experimental progress in determining network parameters, any neural network 
model is underdetermined, because of the sheer complexity of brain tissue and inevitable uncertainties in 
the data. For instance, counting synapses on a large scale presently takes a prohibitive amount of time, and 
no available technique allows determining precise synaptic weights for entire neural populations. Although 
it may be possible to quantify the full connectome of an individual in future, inter-individual variability 
will require modelers to express connectivity as connection probabilities or rules of self-organization if 
they want to learn about general principles of the brain. In modeling studies, parameters are usually 
tuned manually to achieve a satisfactory state of activity, which becomes unfeasible for high-dimensional 
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models due to the size of the parameter space. In particular, it is a priori unclear how parameters of the 
model influence its activity. In consequence, modifications cannot be performed in a targeted fashion, 
and it is difficult to hnd a minimal set of modihcations necessary for aligning a model with experimental 
activity data. 

Overcoming this problem requires a shift of perspective. Instead of regarding the model exclusively 
in a forward manner, generating predicted activity from structure, we in addition consider the system 
in a reverse manner, predicting the structure necessary to explain the observed activity. Our theory, 
starting from a mean-field description, provides a direct link between structure and activity. In contrast 
to simulations, the theory is invertible, which we exploit to identify connections critical for the dynamics 
and to find a minimal set of modifications to the structure yielding a realistic set point of activity. The 
predicted alterations generate hypotheses on brain structure, thus feeding back to anatomical experiments. 

Applying the method to a multi-scale network model of the vision-related areas of macaque cortex, we 
derive targeted modifications of a set of critical connections, bringing the model closer to experimental 
observations of cortical activity. Based on the global properties of the bistable phase space of the model, 
the method refines the model’s construction principles within experimental uncertainties and identifies 
the connections that decisively shape the dynamics. Preliminary results have been presented in abstract 
form (Schmidt et ah, 2015). 

Results 

Global stability in a simple network 

We consider networks with neurons structured into N interconnected populations. A neuron in popula¬ 
tion i receives Kij incoming connections from neurons in population j, each with synaptic efficacy Jij. 
Additionally, each neuron is driven by AText Poisson sources with rate i/ext and synaptic efficacy Jext- 
All neurons in one population have identical parameters, so that we can describe the network activity in 
terms of population-averaged firing rates zzj. 

Our method is applicable if the employed neuron model has a known gain function, either analytically 
or as a function obtained by interpolating numerical results from simulations. In this study, we model 
single cells as leaky integrate-and-fire model (LIF) neurons (see Methods). The possible stationary states 
of these networks are characterized by the firing rates that are equilibria of 

P:=^ = $(«z,A)-iz, (1) 

where s is a pseudo-time. The gain function is known analytically and A indicates its dependence on 
the model parameters {K, J, i^ext, ■ • ■ } (see Methods). 

The input-output relationship $ typically features a non-linearity which, in combination with feedback 
connections, can cause multi-stability in the network. In particular excitatory-excitatory loops cause the 
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system defined by (1) to exhibit multi-stable behavior in the stationary firing rates. A necessary condition 
for the bistability is that the transfer function has an inflection point. The LIF neuron model can have such 
an inflection point, originating from the interplay of its leak term and the threshold. Less realistic neuron 
models, such as the perfect integrate-and-fire model, do not have such an inflection point. To illustrate 
its origin, we first consider the noiseless case (Dayan & Abbott, 2001) without absolute refractoriness 
(ti- = 0). The transfer function initially grows from zero with infinite slope due to the threshold and 
crosses the identity line (Fig. 2C). For larger input the leak term can be neglected and the transfer 
function approaches a linear function with finite slope Ve^Vr (Kriener et ah, 2014), eq. 

11), equivalent to a perfect integrator. This is only possible with a negative curvature at intermediate 
rates, i.e., a reduction in the slope, which makes the transfer function cross the identity line another 
time, causing the bistability. Network-generated noise only affects the low-rate regime where it smears 
out the kink causing the transfer function to grow from zero with positive curvature (see, e.g., (Brunei, 
2000), eq. 22). Importantly, the qualitative picture, i.e., the bistable behavior, is not altered. A finite 
refractory period only has an effect for very high input values where the transfer function saturates at 
I/tj. = 500 spikes/s for the given parameters. 

The basic problem is that there is a trade-off between excitation at the fixed points and their stability. 
In particular, exciting the model to bring a fixed point closer to experimental observations requires a 
method to preserve its stability. We achieve this by controlling the influence of the excitatory-excitatory 
loops on the phase space of the network. 

As an illustration we first study the mechanism using the simple network architecture depicted in 
Fig. 2A: a population of excitatory neurons is coupled to itself with indegree K and is driven by external 
Poisson sources with the same indegree AText = K and rate r'ext • All connections have identical synaptic 
weights J. An increase in the external drive shifts the input-output relationship $(j^, z/ext) of this one¬ 
dimensional system (Fig. 2C) to the left. The bifurcation diagram is shown in Fig. 2E: for low z^ext there 
is only one fixed point with low activity (LA). When increasing iZext, an additional pair of fixed points of 
which one is stable and the other is unstable emerges via a saddle-node bifurcation, leading to a bistable 
system. The second stable fixed point exhibits high firing rates, denoted as the high-activity (HA) state. 
For even higher values of z/ext, the LA state loses stability. 

The equilibria, given by the zeros of the velocity z> in the bistable case, are shown in Fig. 2F. An 
increase of the drive on the one hand increases the firing rate of the LA fixed point (inset) but on the 
other hand reduces its basin of attraction, indicated by the colored bars in the top left corner. For 
illustrative purposes, we extend the problem to two dimensions by splitting the excitatory population 
into two subpopulations of equal size (Fig. 2B), mimicking a loop between excitatory populations in the 
model of the vision-related areas of macaque cortex. The corresponding (symmetric) two-dimensional 
phase space is shown in Fig. 2D. The basin of attraction for the LA fixed point, limited by the separatrix 
(Strogatz, 1994), is reduced with increasing external drive. 

Since we have a bistable system, there must be at least one unstable fixed point on the separatrix at 


4 




Activity constraints on the connectome 


Schuecker, Schmidt et al. 




E Rate V (i/s) 



E z^t'xt (i/s) 



Rate V (i/s) 



0.0 


0.0 0.5 1.0 1.5 


Excitatory rate 2 (10 ^/s) 


Figure 2 Activity flow in an illustrative network example. Left column: Global stability analysis in the 
single-population network. A Illustration of network architecture. C Upper panel: Input-output relationship 
i^ext) for external Poisson drive i/ext = 160 ^ shown in gray. In addition z/ext) for Tr = 0 for the 
noiseless case (blue) and the noisy case (red). The inset shows the gray curve over a larger input range. Lower 
panel: $(zy, z^ext) for different rates of the external Poisson drive z^ext = [150,160,170] ^ from black to light 
gray. Intersections with the identity line (dashed) mark fixed points of the system, which are shown in E as a 
function of z/gxt- F Flux z> in the bistable case for $(z/ext = 1601, A") in black, <I’(z/'xt = 16ll, A") in blue, 
and modified system ^‘(I'ext = 16ll,Ar') in red. Intersections with zero (dashed line) mark fixed points. The 
inset shows an enlargement close to the LA fixed point. Horizontal bars at top of figure denote the size of 
the basin of attraction for each of the three settings. Right column: Global stability analysis in the network 
of two mutually coupled excitatory populations. B Illustration of network architecture. D Flow field and 
nullclines (dashed curves) for ^•(r'ext = 160-1, A') and separatrices (solid lines), LA fixed point (rectangle), 
HA fixed point (cross) and unstable fixed points (circles) for $(z/ext = 160-1, AT) in black, ^(t'ext = 161-1, AT) 
in blue, and #(z/ext = 16ll, A'') in red. The red separatrix and the red unstable fixed point coincide with the 
black ones. G Enlargement of D close to the LA fixed points. Flow field of original system shown in black, of 
modified system in red. 
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the intersection of the nullclines, i.e., the subspace for which the velocity Vi in direction i vanishes. We 
use the unstable fixed point to preserve the basin of attraction when the external drive t'ext is increased. 
For this purpose, we modify the connectivity K —>■ K' to reverse the shift of the unstable hxed point due 
to the parameter change r'ext —t r'ext (see Methods for a detailed derivation). Since the separatrix follows 
the unstable fixed point, this approximately restores the original basin of attraction. 

The resulting velocity of the system ($ defines the system (1)) is shown in Fig. 2F. The 

firing rate in the LA fixed point is increased as desired (inset), and the unstable fixed point coincides with 
that obtained in the original system <l>(r'ext, K). This pattern of fixed points is also indicated by the zero 
vectors of the velocity field v (Fig. 2D). The separatrix follows the unstable fixed point, and the basin of 
attraction in the system is restored to that in the original system ^(r'ext, AT). Fig. 2G shows 

the behavior of the LA fixed point in more detail. The modihcation of K does not noticeably change the 
location of the LA fixed point. In conclusion, the method allows us to increase the firing rates in the LA 
fixed point without modifying its basin of attraction. 

The purely excitatory network is the simplest model to explain a phase space configuration with a LA 
and a HA fixed point. Inhibitory feedback is not necessary for this bistability, but it would certainly alter 
the input-output relationship. For example the classical excitatory-inhibitory network (Brunei, 2000) 
in the balanced regime has an input-output relationship with a negative slope and thus only one fixed 
point exists. However, if a pair of such balanced E-I networks is coupled by sufficiently strong mutually 
excitatory connections, these connections cause a bistability in a similar manner. Thus the mechanism 
shown in the purely excitatory network can also lead to the emergence of a HA attractor in more complex 
networks with inhibition. 


Bistability in the multi-scale network model 

We investigate a multi-scale network model of cortical areas to understand the structural features essential 
for a realistic state of baseline activity. The model extends and adapts the microcircuit model presented 
in (Potjans & Diesmann, 2014), which covers Imm^ surface area of early sensory cortex (Fig. 3A), to all 
vision-related areas of macaque cortex (Fig. 3B). Based on the microcircuit model, an area is composed 
of 4 layers (2/3, 4, 5, and 6) each having an excitatory (E) and an inhibitory (I) neural population, except 
parahippocampal area TH, which consists of only 3 layers (2/3, 5, and 6). A detailed description of the 
data integration is given in (Schmidt et ah, 2016). 

Simulations of the model (Fig. 3C) reveal that, though realistic levels of activity can be achieved for 
populations in layers 2/3 and 4, populations 5E and 6E of the majority of areas show vanishingly low 
or zero activity in contrast to empirical data (Swadlow, 1988; de Kock & Sakmann, 2009). Inputs from 
subcortical and non-visual cortical areas are modeled as Poissonian spike trains, whose rate t'ext is a free, 
global parameter. To elevate the firing rates in the excitatory populations in layers 5 and 6, we enhance 
the external Poisson drive onto these populations parametrized by n (see Methods). However, already 
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Figure 3 Bistability of the model. A Sketch of the microcircuit model serving as a prototype for the areas of 
the multi-area model (figure and legend adapted from figure 1 of Potjans and Diesmann (Potjans &l Diesmann, 
2014), with permission). B Sketch of the most common laminar patterns of cortico-cortical connectivity of 
the multi-area model. C Population-averaged firing rates encoded in color for a spiking network simulation of 
the multi-area model with low external drive {k = 1.0). D As C but for increased external drive {k = 1.125). 
The color bar refers to both panels. Areas are ordered according to their architectural type along the horizontal 
axis from VI (type 8) to TH (type 2) and populations are stacked along the vertical axis. The two missing 
populations 4E and 41 of area TH are marked in black and firing rates < 10“^ Hz in gray. E Histogram 
of population-averaged firing rates shown in C for excitatory (blue) and inhibitory (red) populations. The 
horizontal axis is split into linear- (left) and log-scaled (right) ranges. F as E corresponding to state shown in 
D. 
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a perturbation of a few percent leads to a state with unrealistically high rates (Fig. 3D), caused by the 
reduced basin of attraction of the low-activity state similar to Fig. 2D. Our aim is to improve the working 
point of the model such that all populations exhibit spiking activity > 0.05 spikes/s while preventing 
the model from entering a state with unrealistically high rates of > 30 spikes/s (figure 13 of (Swadlow, 
1988), (de Kock & Sakmann, 2009)). The employed technique exposes the mechanism giving rise to the 
observed instability and identifies the circuitry responsible for this dynamical feature. 

Targeted modifications preserve global stability 

We apply the procedure derived in Methods and find targeted modifications to the connectivity K that 
preserve the global stability of the low-activity fixed point for increased values of the external drive k. 

In the following we choose the inactive state ^^(O) = (0,... ,0)"^ as the initial condition. The exact 
choice is not essential since we are only interested in the fixed points of the system. Fig. 4B shows the 
integration of (1) over pseudo-time s for different levels of the external drive to populations 5E and 6E 
parametrized by k. Eor low values of k, the integration converges to the LA fixed point shown in Eig. 4D, 
and is in agreement with the activity emerging in the simulation (Fig. 3C). For increased values of k, the 
system settles in the HA fixed point (Fig. 4E), again in agreement with the simulation. The population- 
specific firing rates in the HA state found in the mean-field predictions (Fig. 4E) are also close to those in 
the simulation (Fig. 3D), but minor deviations occur due to the violation of the assumptions made in the 
diffusion approximation. In particular, at these pathologically high rates, the neurons fire regularly and 
close to the reciprocal of their refractory period, while in the mean-field theory we assume Poisson spike 
statistics. Still, the mean-field theory predicts the bistability found in the simulation. Since the theory 
yields reliable predictions in both stable fixed points, we assume that also the location of the unstable 
fixed point in between these two extremes is accurately predicted by the theory. 

To control the separatrix we need to find the unstable fixed point of the system. This is nontrivial 
since the numerical integration of (1) for finding equilibria by construction only converges to stable fixed 
points. If the unstable fixed point has only one repelling direction (Fig. 5A), it constitutes a stable 
attractor on the N dimensional separatrix. The separatrix is a stable manifold (Strogatz, 1994), and 
therefore a trajectory originating in its vicinity but not near an unstable fixed point initially stays in the 
neighborhood. If an initial condition just outside the separatrix is close to the basin of attraction of a 
particular unstable fixed point, the trajectory initially approaches the latter. Close to the fixed point 
the velocity is small. Ultimately trajectories diverge from the separatrix in the fixed point’s unstable 
direction, as illustrated in Fig. 4A. In conclusion, we expect a local minimum in the velocity along the 
trajectories close to the unstable fixed point. To estimate the location of the unstable hxed point in this 
manner, we need to find initial conditions close to the separatrix. Naively, we would just fix the value of k 
and vary the initial condition. However, due to the high dimensionality of our system this is not feasible 
in practice. Instead, we vary k, for a fixed initial condition. Fig. 4B shows the firing rate averaged across 
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Figure 4 Application of the mean-field theory to the multi-area model. A Trajectories of (1) starting 
inside the separatrix converge to an unstable fixed point. Trajectories starting close to the separatrix are 
initially attracted by the unstable fixed point but then repelled in the fixed point's unstable direction and 
finally converge to a stable fixed point. B Firing rate averaged across populations over time. Integration of 
(1) leads to convergence to either the low-activity (LA) or the high-activity (HA) attractor for different choices 
of the external input factor k, with k = 1 the original level of external drive. We show eight curves with n 
varying from 1.0 to 1.014 in steps of 0.002 and two additional curves for k = 1.007662217, 1.007662218. The 
curves for the largest factor (k = 1.007662217) that still leads to the LA state and for the smallest factor 
(k = 1.007662218) that leads to the HA state are marked in blue and red, respectively. The four curves with 
K < 1.006 coincide with the blue curve. C Euclidean norm of the velocity vector in the integration of (1) for 
the different choices of k. The vertical dashed line indicates the time Sc of the last local minimum in the blue 
curve. D Stationary firing rate in the different areas and layers of the model in a low-activity state for k = 1.0 
as predicted by the mean-field theory (same display as in Fig. 3). E As D, but showing the high-activity state 
for K = 1.125. 
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Figure 5 Eigenspectrum analysis of network stability. A Eigenvalue spectrum of the effective connectivity 
matrix M for the first (blue squares) and second (red dots) iteration. The dashed vertical line marks the 
edge of stability at a real part of 1. B Contribution rji ((13)) of an individual eigenprojection I to the shift of 
the unstable fixed point versus the relative change in indegrees associated with I for the first iteration. The 
data point corresponding to is marked in red. The inset shows the relative angles between Sv>* and the 
eigenvectors u^. The red line corresponds to the critical eigendirection. C Entries of the eigenvector 
associated with in the populations of the model. The affected areas are 46 and FEE. D Same as C for 
the second iteration. 



Figure 6 Unstable fixed points in subsequent iterations. Population firing rates at the unstable fixed 
point as predicted by the mean-field theory encoded in color for iterations 1 (A) and 2 (B). Same display as 
in Fig. 3. 


populations for two trajectories starting close to the separatrix, where the first one converges to the LA 
fixed point and the second one to the HA state. The trajectories diverge near the unstable fixed point 
and thus we define the last local minimum of the Euclidean norm of the velocity vector as the critical 
time Sc at which we assume the system to be close to the unstable fixed point (Fig. 4C). We find four 
relevant and distinct unstable fixed points, of which two are shown in Fig. 6. 

To counteract the shift of the separatrix caused by the increase in k, we follow the procedure described 
in Methods. We subject the modifications of connectivity to the additional following constraints. In 
line with the anatomical literature, we do not allow for changes of the connectivity that would lead to 
cortico-cortical connections originating in the granular layer 4 (Felleman & Afan Essen, 1991), and we 
also disallow inhibitory cortico-cortical connections, as the vast majority of long-range connections are 
known to be excitatory (Salin & Bullier, 1995; Tomioka & Rockland, 2007). In addition, we naturally 
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restrict indegrees to positive values. We find that four iterations (numbered by index j) corresponding 
to the four distinct unstable fixed points suffice to preserve the basin of attraction of the LA state with 
respect to an increase of the external drive up to k = 1.15. In the following we concentrate on iterations 
1 and 2, where the second one is also representative for iterations 3 and 4, which are qualitatively alike. 
To derive the required modifications of the indegree matrix, we decompose K into its N eigenmodes 
and quantify the contribution of each eigenmode to the shift of the unstable fixed point (see Methods). 
This allows us to identify the most effective eigendirection: in each iteration j there is exactly one 
unstable eigendirection with an eigenvalue Re > 1 (Fig. 5A). The associated critical eigenvector is 

approximately anti-parallel to the shift of the fixed point, Su* (inset of Fig. 5B), and of similar length. 
The critical eigendirection (red dot in Fig. 5B) constitutes the most effective modification, giving the 
largest contribution to the desired shift while requiring only a small change of 2.3% in average total 
indegrees. In the chosen space of eigenmodes, the modifications are minimal in the sense that only this 
most effective eigenmode is changed. 

The associated eigenvector predominately points into the direction of populations 4E and 5E of 
areas FEF and 46 (Fig. 5C), while Uc has large entries in the 5E populations of two areas (Fig. 5D). The 
high rates of these populations at the unstable fixed points (cf. Fig. 6A,B with Fig. 5C,D) reflect that 
the instability is caused by increased rates in excitatory populations, particularly in population 5E. Each 
iteration shifts the transition to the HA state (the value of k for which the separatrix crosses the initial 
condition) to higher values of k and increases the attainable rates of populations 5E and 6E in the LA 
state (Eig. 7A). After all four iterations, the average total indegrees (summed over source populations) 
of the system are changed by 11.3%. The first iteration mainly affects connections within and between 
areas 46 and FEE (Eig. 7B). In particular, the excitatory loops between the two areas are reduced in 
strength, especially those involving layer 5 (Fig. 7C). We thus identify two areas forming a critical loop. 
In the remaining iterations, the changes are spread across areas and especially connections originating 
in layer 5 are weakened (Fig. 7D). In conclusion, the method identifies critical structures in the model 
both on the level of areas and on the level of layers and populations, and leads to a small but specific 
structural change of the model. 

Analysis of the modifications 

In the following we analyze the modifications of the connectivity with respect to the internal and inter¬ 
area connections in detail. The intrinsic circuits of the areas are modified in different directions, as 
shown for two exemplary areas V4 and CITv in Fig. 8A. Despite this heterogeneity, significant changes 
affect mostly excitatory-excitatory connections (Fig. 8A, bottom panel) with connections from popula¬ 
tion 5E experiencing the most significant changes (top panel of Fig. 8A). In fact, the anatomical data 
(Binzegger et ah, 2004) underlying the microcircuit model (Potjans & Diesmann, 2014) contain only 
two reconstructed excitatory cells from layer 5, but considerably more for other cell types, indicating a 
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Figure 7 Altered phase space and modified connections. A Firing rates averaged across populations 5E 
and 6E and across areas for different stages from the original model (black) to iteration 4 (light gray) as a 
function of K, predicted by the mean-field theory. B Relative changes in the indegree AT^b between 

areas A^B in the first iteration. C Layer-specific relative changes 5Kin/ 'Y^^Kin in the connections within 
and between areas FEF and 46, for the first iteration. Populations are ordered from 2/3E (left) to 61 (right) on 
the horizontal axis and from 61 (bottom) to 2/3E (top) on the vertical axis as in panel D. D Relative changes 
in population-specific indegrees summed over target populations, combined for iterations 

two, three and four. 


higher uncertainty for layer 5 connections. Fig. 8B shows the correlation between intrinsic connectivity 
changes for all pairs of areas, with areas ordered according to a hierarchical clustering using a farthest 
point algorithm (Voorhees, 1986) on the correlation matrix. We find four clusters each indicating a group 
of areas which undergo changes with similar patterns. The groups are displayed in different colors in 
the histogram in Fig. 8B. The areas of the model are categorized into architectural types based on cell 
densities and laminar thicknesses (see Methods). Areas with architectural type 4, 5 and 6 are distributed 
over several clusters. We can interpret this as a differentiation of these types into further subtypes. 
The resulting changes of the intra-areal connectivity are small (Fig. 8A), but still significant for network 
stability. 

Connections between areas can be characterized by their FLN and SLN (see Methods). The FLN 
reflects the overall strength of an inter-areal connection and is only weakly affected across connections 
(Fig. 8C), with a correlation between original and modified logarithms of FLN of rpearson = 0.79. Sig¬ 
nificant variations in the FLN occur mostly for very weak connections that are likely to have substantial 
relative uncertainties in the experimental data. The two overlapping red dots in Fig. 8C represent the 
connections between areas 46 and FEF, which are modified in the first iteration (Fig. 7C). The SLN 
determines the laminar pattern of the location of source neurons for cortico-cortical connections. Over¬ 
all, data are available for 24% of the inter-areal connections in the parcellation of Felleman & van Essen 
(Felleman & Van Essen, 1991), while the SLN for the rest are derived from the sigmoidal law. The 
majority of connections undergo small changes in their laminar source pattern (Fig. 8D) and connections 
with large modifications (1(5577^1 > 0.5) are weak (average FLN = 6-10“^ compared to FLN = 10“^ in 
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SSLN 

Figure 8 Analysis of changes in connectivity. A Top panel: relative changes in population-specific intrinsic 
indegrees summed over target populations and averaged across areas, ^^aaI Si ^^'aa)a- Bottom panel: 
changes in the indegrees within and between exemplary areas V4 and CITv relative to the total indegrees of the 
target populations, i.e., Populations ordered as in Fig. 7C. B Pearson correlation coefficient 

of the changes of the internal indegrees between all pairs of the 32 areas. Areas ordered according to 

hierarchical clustering using a farthest point algorithm (Voorhees, 1986). The heights of the bars on top of 
the matrix indicate the architectural types of the areas (types 1 and 3 do not appear in the model) with color 
representing the respective clusters. C FLN of the modified connectivity after 4 iterations versus the original 
FLN of the model. Only FLN > 10“^ are shown for a better overview. The overlapping red dots represent 
the connections between areas 46 and FEF. Unity line shown in gray. D Histogram of the cumulative changes 
in SLN over all four iterations {SSLN = SLN^'^^ - SLN<^°'>). 
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the model as a whole). Because weak connections are represented by low counts of labeled neurons, they 
have a relatively large uncertainty in their laminar patterns, justifying larger adjustments. Spearman’s 
rank correlation between the SLN of the original model that were directly taken from experiments and 
the logarithmic ratios of cell densities is p = —0.63 (p = 3 • 10“^^, p-value of a two-sided test for un¬ 
correlated data). For the modified model, we take the SLN of all connections into account and obtain 
p = —0.40 (p = 6-10“^°), indicating a reduced, but still significant, monotonic dependence between SLN 
and the logarithmic ratios of cell densities. 

To judge the size of the modification to the connectivity, we compare it to the variability of measured 
cortico-cortical connection densities (Markov et al., 2014a). We quantify the latter as the average inter¬ 
individual standard deviation of the logarithmic FLN, i.e., a = ^ (log FLN — log FLN ^"^, where the 
overbar ” denotes the average over injections and (•) the average over connections. This variability equals 
2.17 while the average modification of the logarithmic FLN is 1.34. The main experimental connection 
probabilities used to construct the intra-areal connectivity of the model have an average relative standard 
deviation of 30% across electrophysiological experiments (cf. Table 1 of (Potjans & Diesmann, 2014)) 
while the intra-areal connection probabilities of the model are modified by 9% on average. The authors of 
(Scannell et ah, 2000) report even greater variability in their review on cortico-cortical and thalamocortical 
connectivity. These considerations show that on average, the changes applied to the connectivity are well 
within the uncertainties of the data. Overall, 35 out of 603 connections were removed from the network. 
In the CoCoMac database, 83 % of these are indicated by only a single tracer injection, while the overall 
proportion of connections measured by a single injection is 59%. 

For the modified connectivity and n = 1.125, which we choose to avoid being too close to the transition 
(Fig. 7A), the theory predicts average rates in populations 5E and 6E of 1.3 and 0.18 spikes/s, which is 
closer to experimentally observed rates compared to the original model. Furthermore we find that the 
modified connectivity allows us to decrease the inhibition in the network to g = —11. Simulating the full 
spiking network model then results in reasonable rates across populations and areas (Fig. 9B, D). The 
average rates in populations 5E and 6E are increased compared to a simulation of the original model 
from 0.09 and 2 • 10“® spikes/s to 3.0 and 0.4 spikes/s, respectively. All populations exhibit firing rates 
within a reasonable range of 0.05 to 30 spikes/s (Fig. 9D), as opposed to the original state in which a 
considerable fraction of excitatory neurons is silent (Fig. 3E). The theoretical prediction is in excellent 
agreement with the rates obtained in the simulation (Fig. 9A, C). Small discrepancies are caused by 
violations of the employed assumptions, i.e., Poissonian spiking statistics (Fourcaud & Brunei, 2002). 
Differences between theory and simulation are small, and negligible for the central aim of the study: the 
integration of activity constraints into the data-driven construction of multi-scale neuronal networks. 
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Figure 9 Improved low-activity fixed point of the model. Population-averaged firing rates for k = 1.125 
encoded in color (A) predicted by the analytical theory and (B) obtained from the full simulation of the spiking 
network. Same display as in Fig. 3. C Analytical versus simulated firing rates (black dots) and identity line 
(gray). D Histogram of population-averaged simulated firing rates. Same display as in Fig. 3. 


Discussion 


This study investigates the link between experimentally measured structural connectivity and neuronal 
activity in a multi-scale spiking network model of the vision-related areas of macaque cortex (Schmidt 
et al., 2016). We devise a theoretical method that systematically combines anatomical connectivity with 
physiological activity constraints. Already weak constraints, demanding the activity to neither vanish nor 
be pathologically high, yield a set of specific but small structural modihcations necessary to increase the 
model’s excitation to a realistic level. We do not fit the model parameters to experimental activity data 
in the sense of minimizing an error function, since the sparse relevant experimental data do not allow for 
defining such a function in a meaningful way. Nevertheless, we considerably improve the network activity 
with a change in only a small set of parameters. The procedure constrains the experimentally obtained 
connectivity maps to a realization that is compatible with physiological experiments. This establishes a 
path from experimentally observed activity to specific hypotheses about the anatomy and demands for 
further experiments. 

Connections are modified both inside and between cortical areas, on average well within the known 
uncertainties of the underlying data. The model areas are based on the early sensory cortex model pre¬ 
sented in (Potjans & Diesmann, 2014). This circuit is adapted to individual areas by taking into account 
neuronal densities and laminar thicknesses. The model definition renders areas with equal architectural 
type similar in their internal connectivity, a drastic but inevitable simplihcation due to the lack of more 
detailed experimental data. The proposed method softens this assumption: small adaptations of internal 
connectivity distinguish the architectural types into further subtypes. These modifications are signif¬ 
icant for the global stability of the network. Thus, our approach enables purely anatomy-based area 
categorizations to be refined with dynamical information. 
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Connections between areas are changed in terms of total strength and laminar patterns. Overall, the 
changes are small, but significant for specific connections. The loop formed by areas 46 and FEF is critical 
to the global stability of the network. Both areas have been investigated in (Markov et ah, 2014a), albeit 
in a different parcellation scheme than the scheme used here (Felleman & Van Essen, 1991). Our method 
suggests a weaker coupling of these two areas than found in the anatomical data set. Uncertainties, 
partly due to the mapping between parcellations, leave room for this interpretation. Areas 46 and FEF 
belong to prefrontal cortex and are multimodal, indicating that the influence of other parts of cortex 
could stabilize them, a mechanism outside the scope of the present model of vision-related areas. Both 
explanations can be tested either in an experimental study or with an extended model. 

A few weak connections undergo large changes in their laminar patterns. With the present activity 
constraints, the method hereby weakens hierarchical relations in the model structure, indicating that pre¬ 
serving these relations requires additional dynamical constraints such as layer-specific coherence between 
areas (?Bastos et ah, 2015). Conversely, it is possible to achieve satisfactory population rates with a less 
pronounced hierarchical structure. 

Our analysis reveals that layer 5 excitatory cells play a critical role in the model’s dynamics, in 
line with the observed ongoing activity in mammalian neocortex (Sanchez-Vives & McCormick, 2000; 
Beltramo et ah, 2013). This critical role is often attributed to single-neuron properties, with a subset of 
layer 5 neurons displaying pacemaker activity (Le Bon-Jego & Yuste, 2007; Lorincz et ah, 2015; Neske 
et ah, 2015). In addition, we here find that the network architecture itself already explains the strong 
impact of layer 5 on the phase space of the network, suggesting that single-neuron properties and network 
structure jointly enable layer 5 to exert its dominant influence. 

The cortico-cortical connectivity of the model is compiled from the extensive dataset of (Markov 
et ah, 2014a) combined with the CoCoMac database (Stephan et ah, 2001; Bakker et ah, 2012), which 
collects data from hundreds of tracing studies. One could consider alternative methods for combining 
this information into one connectivity graph, for instance taking into account how consistently a given 
connection is reported across studies (Schmitt et ah, 2014), and compare different methods by analyzing 
the resulting network dynamics. The presented mean-held theory could then be used to estimate the 
bring rates of each network instance without performing time-consuming simulations. However, here we 
hrst choose the integrative approach that accumulates all available experimental evidence into a single 
model and afterwards modify the resulting connectivity with our analytical procedure, thereby effectively 
discarding uncertain connections. 

We restrict this study to networks of leaky integrate-and-hre model neurons, consistent with the key 
concept of the models we consider: individual cells are modeled in a simple manner to expose the impact 
of structural connectivity on the network dynamics. Moreover, the current-based leaky integrate-and-hre 
neuron can reproduce in-vivo like activity (Rauch et al., 2003; Jolivet et ah, 2006) and is analytically 
tractable, which enables the identihcation of mechanisms underlying specihc network effects. More com¬ 
plex neuron models can be incorporated into the method by replacing the gain function of each neuronal 
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population with an analytical expression or an interpolated function obtained from spiking single-neuron 
simulations. For example, one could use a conductance-based point-neuron model for which the network 
dynamics can be described by population rate models (Shriki et ah, 2003), featuring a non-monotonic 
gain function: the gain is reduced if excitatory and inhibitory inputs are increased in a balanced manner 
(Kuhn et ah, 2004). Generally, this renders a system more stable. However, the bistability considered 
in our work is caused by excitatory inputs. Since conductance-based models also have a monotonically 
increasing gain function in dependence of the excitatory conductance alone, we expect the bistability to 
occur for such models as well. 

Importantly, our method only requires a description of the system’s fixed points and their dependence 
on model parameters. The employed theory uses the diffusion approximation to derive a self-consistency 
equation for the stationary population rates valid for high indegrees and small synaptic efficacies. These 
requirements are fulfilled in the multi-area model and therefore the theoretical prediction agrees with the 
stationary activity of the simulation. Moreover, the theory predicts the bistability of the model, which is 
non-trivial, as the mean-field assumption of Poisson statistics of the activity is generally violated in the 
high-activity state. Nevertheless, since the activity in this state is mostly driven by strong mean inputs 
and the theory converges to the noiseless solution in the mean-driven limit, its predictions still provide 
viable approximations. The firing rates in the unstable fixed points are predominantly low, while the 
exceptions with very high rates are again mean-dominated. Presumably this explains the accuracy by 
which the theory predicts the locations of these fixed points and the resulting global stability properties. 

Since the high-activity attractor of the model under consideration is unrealistic, we aim to prevent a 
transition to this state. However, in the high-dimensional system it is not a priori clear in which direction 
the separatrix has to be shifted to ensure stable dynamics. We therefore choose a pragmatic approach and 
shift the separatrix back to its initial location, inverting the shift which reduced the global stability of the 
low-activity fixed point. We achieve this to a good approximation by preserving the location of unstable 
fixed points on the separatrix. To this end, we use a linearization around these locations and an eigenmode 
decomposition to identify the set of connections to adjust. In the multi-area model, this linearization is 
justified because the system operates close to an instability so that only minor modihcations are required. 
The method can be generalized to larger modihcations by changing parameters iteratively in small steps. 

Biological networks have various stabilization mechanisms not considered here, which render them 
less critical. For instance, during growth, homeostatic mechanisms guide the system toward the right 
structure. Furthermore, short-term synaptic plasticity (reviewed in (Morrison et ah, 2008)), homeostatic 
synaptic scaling (Turrigiano et ah, 1998) and spike-frequency adaptation (e.g. summarized in (Benda & 
Herz, 2003)) may prevent the system from entering the high-activity state. However, introducing these 
self-organizing mechanisms increases model complexity, causing a more intricate relation between struc¬ 
ture and activity. Therefore, we start from a mean-held description on the level of neuronal populations, 
ignoring details of synaptic dynamics. Mild constraints on the activity lead to a network structure within 
the anatomical range of parameters. This network yields globally stable activity, suggesting that addi- 
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tional stabilization mechanisms are not required to achieve this. Nonetheless, they can potentially render 
the network more robust against external stimulation. 

The observed inter-individual variability may reflect that mechanisms of self-organization and home¬ 
ostasis find structurally different implementations of the same function. Thus, in studies across individuals 
we cannot expect that progress in experimental technology narrows down the variability of parameters in¬ 
definitely. The combined ranges of parameters rather specify the solution space, and our method provides 
a way to find a particular solution. 

In principle, the method applies to any model parameter. It would be possible to modify, for example, 
synaptic weights. Since experimental data on synaptic weights are sparse, this is another natural choice. 
Moreover, such an analysis may provide hints about suitable synaptic plasticity rules that dynamically 
stabilize the model. The method can be applied to networks with more complex sets of attractors 
compared to the bistable case considered here. Though in high-dimensional systems such as our multi¬ 
area model, a larger number of attractors would make it more challenging to find all relevant unstable 
fixed points, the underlying idea of preserving the location of a separatrix is general. In contrast to the 
model considered here, transitions between fixed points can have a functional meaning in certain neuronal 
networks with multiple attractors. The specific location of the separatrix is then functionally relevant. 
Our method exposes the sensitivity of the location of the separatrix to certain model parameters and 
allows controlling its location in a specific manner. 

In this work, we analyzed the global stability properties of the neuronal network on the population 
level. In contrast, Ostojic (Ostojic, 2014) performs a local stability analysis on the level of single neu¬ 
rons of an initially stable fixed point in a system with only one attractor. The author investigates the 
point at which the real parts of the eigenvalues of the Jacobian matrix evaluated at this fixed point 
become positive, i.e., the fixed point turns spectrally unstable and the system undergoes a transition 
to a heterogeneous asynchronous state. Analyzing the spectral stability on the single-neuron level does 
not reveal the global stability properties required in the current work: While a local stability analysis 
only considers infinitesimal perturbations, studying the basin of attraction gives information about the 
size of fluctuations against which the fixed point is stable. We expect both attractors to be spectrally 
stable because they do not show strong rate fluctuations and the mean-field theory predicts the activity 
accurately in both cases. Nevertheless, a heterogeneous state could occur if the synaptic weights were 
increased or the external drive was stronger. However, (Goedeke et ah, 2016) show that the transition 
in stochastic systems that quantitatively resemble the spiking network (Grytskyy et ah, 2013), does not 
coincide with the loss of spectral stability. 

One striking feature of the heterogeneous state is bursty spiking behavior of individual cells. Bursty 
spiking is also observed in the multi-area model for increased synaptic weights of inter-area connections 
(Schmidt et ah, 2016). The fixed points cannot be accurately described in this case because the fluctua¬ 
tions need to be taken into account (Mastroguiseppe & Ostojic, 2016). Simulations show (figures 4 and 
5 of (Schmidt et ah, 2016)) that the modifications obtained in this study are still able to prevent the 
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system from a transition to a HA attractor also in the presence of bursting neurons. This indicates that 
the phase space does not change qualitatively and our results are robust against such bursting behavior. 

Experimental data on stationary activity in cortex are sparse. We therefore restrict ourselves to 
fundamental constraints and increase the drive to the model in an area-unspecific way to fulfill them. 
The resulting heterogeneity of the firing rates across areas is thus not imposed by the method, but rather 
arises from the connectivity that remains strongly informed by anatomical data. Alternatively, one could 
predefine a desired state and investigate the parameter changes necessary to achieve it. 

The presented analytical method that combines anatomy and activity data into a consistent model 
is restricted to stationary firing rates. In future studies, also higher-order statistical measures of activity 
can be used as constraints. Resting-state fMRI, for example, provides information on the functional 
connectivity between areas as a second-order measure. When combined with analytical predictions of 
functional connectivity, the method may shed light on the anatomical connection patterns underlying 
inter-area communication. 


Methods 


In this study, we model single cells as leaky integrate-and-fire model neurons with exponentially decaying 
postsynaptic currents. Table 1 specifies the model parameters. In the diffusion approximation, the 

Table 1 Specification of the neuron and synapse parameters. 


Synapse parameters 

Name 

Value 

Description 

J±SJ 

9 

87.8 ± 8.8 pA 
-16 (Fig. 3) 

-11 (Fig. 9) 

excitatory synaptic strength 
relative inhibitory synaptic strength 

de ± 6de 

1.5 ±0.75 ms 

local excitatory transmission delay 

di ± Sdi 

0.75 ± 0.375 ms 

local inhibitory transmission delay 

d± 6d 

d = s/vt ± \s/vt 

inter-areal transmission delay, with s the 
distance between areas 

Vt 

3.5 m/s 

transmission speed 

Neuron model 

Name 

Value 

Description 

Tm 

10 ms 

membrane time constant 

Tv 

2 ms 

absolute refractory period 

Ts 

0.5 ms 

postsynaptic current time constant 

Cvn 

250 pF 

membrane capacity 

Vv 

-65 mV 

reset potential 

9 

-50 mV 

fixed firing threshold 

Ei^ 

-65 mV 

leak potential 


dynamics of the membrane potential V and synaptic current Is is (Fourcaud & Brunei, 2002) 


T-m-TT = -V + Is(t) 
at 

Ts-yy = -4 + M + (2) 
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where Tm is the membrane time constant and Tg the synaptic time constant, respectively. The membrane 
resistance has been absorbed into the definition of the current. The input spike trains are 

approximated by a white noise current with fluctuations oc and mean value /i. Here ^ is a centered 
Gaussian white process satisfying (^(t)) = 0 and — Whenever the membrane potential 

V crosses the threshold 0, the neuron emits a spike and V is reset to the potential Vr, where it is clamped 
during All neurons in one population have identical parameters, so that we can describe the network 
activity in terms of population-averaged firing rates that depend on population-dependent input fii, (Ji 
determined by the connectivity. Using the Fokker-Planck formalism, the stationary firing rates for each 
population i are given by (Fourcaud & Brunei, 2002) 


1 

6 — (A) 1 _ / Ts 



= Tr+Tmy/TT (I -|-erf(x)) dx 


Vi 

cr^ (A) ^ ’ V "^m 



=: l/$,(iz,A) 

(3) 


— ^ ^ ^ ‘i-j ^m-^ext'^ext^ext 

(4) 

<^Ua) 


(5) 




which is correct up to linear order in and where 7 = \Q{l/2)\/^/2, with ^ denoting the Rie- 

mann zeta function (Abramowitz & Stegun, 1974). Here, A is chosen from the set of model param¬ 
eters {K^J. z/ext,---}. If A is a matrix, we vectorize it by concatenating its rows and indicate this 
by lower case, i.e., a — (ooo, ooij ■ • ■; goat, oiq, . ■., oiat, Qno • ■ •; o,nn) = vec(A"^)'^ following (Magnus & 
Neudecker, 1995). If the chosen parameter is a scalar we denote it with a. 

We find the fixed points of (3) by solving the first-order differential equation (I) (Wong & Wang, 
2006) 


i> := — = $(zz, A) - iz, 

using the classical fourth-order Runge-Kutta method (RK4) with step size h = O.OI, where s denotes a 
dimensionless pseudo-time. The same approach can be used to solve the activity on a single neuron level 
(Sadeh & Rotter, 2015). Note that (1) does not reflect the real time evolution of the population rates, 
but rather is a mathematical method to obtain the system’s fixed points. In contrast to (Wong & Wang, 
2006) we do not only search for stable fixed points, but also use ( 1 ) to obtain unstable attractors (cf. 
Results), an idea originating from the study of simple attractor networks ((Amit & Brunei, 1997b) esp. 
their figure 2 and eq. 12 ). 

In a bistable situation, the initial condition of (1) determines which fixed point the system settles 
in. However, studying the behavior for a particular initial condition is of minor interest, since the actual 
spiking network is a stochastic system which fluctuates around the fixed points of the deterministic 
system defined by (I). Even if we knew that (1) would relax to the LA fixed point for one particular 
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initial condition, this would not necessarily imply that this state is indefinitely stable. Global stability is 
determined by the size of the basin of attraction of the LA fixed point. 

In the following, we derive the equations leading us to targeted modifications of a parameter b neces¬ 
sary to compensate for the changes in the global stability induced by the change of another parameter a. 
To this end, we study the behavior of the fixed points with respect to an infinitesimal change Sa — a' — a 
in the chosen model parameter. Let v'*{a) and v*{a') be the corresponding locations of the fixed points 
and dv* — u*{a') — their separation. We can then expand v*{a') into a Taylor series up to first 

order in 5a and obtain 


with 


u*{a') = u*{a) + ^a5a 

5u* = Ao^a, 


^a^ij — 


dvj 

dttj 

dttj 

d^i doj 
Si 


d^i 1 dgf 
d(7i 2ai daj ’ 


( 6 ) 


(7) 


where we notice that Si and Ti only depend on the target population i. We accordingly define two diagonal 
matrices S and T with Sa = Si and Ta = Ti. We further define the connectivity matrix W = If ® J, 
where @ denotes element-wise multiplication, also called the Hadamard product (see (Cichocki et ah, 
2009), for a consistent set of symbols for operations on matrices). The derivatives with respect to aj have 
the compact expressions 


d^ij _ dUi ^ dni dvn 
duj duj ^ dvn daj 

” =Wi„ 

with the Jacobian (D^f)^^- := of some vector-valued function f and 

d. ^ 

= (Do,<T -L Tm ^^ {KinJin)^a,nj , 

where we use the Hadamard product again to define the matrix W 2 := K ® J @ J. Inserting the total 
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derivatives into (7), we derive the final expression for Aa, reading 


Aa = S [Da/X + TmWA^] + T [DaCr^ + TmW 2 A^ 


Aa = 


t-T^{SW + TW 2 ) 


-.M 


[5Da/x + TT>a.cr 


— Aa 


( 8 ) 


where we use 1 for the identity matrix and define the effective connectivity matrix M and the matrix 
Aa- The latter has dimensionality N x P, where P is the dimension of a (for example, P = N'^ for 
a = k the vector of indegrees). With the aid of (8), evaluating (6) at the unstable fixed point predicts 
the shift of the separatrix (Fig. 2D) to linear order. We now consider an additional parameter b which is 
modified to counteract the shift of the unstable fixed point caused by the change in parameter a, i.e., 


AaSa 

AaSa 


-AtSb 


-AbSb, 


(9) 


where we used that the inverse of 1 — M appears on both sides of the equation and hence drops out. 
Note that the tuple (a, b) may represent any combination of model parameters, for example external 
input, indegrees, synaptic weights, etc. For a particular choice of (a, b) we solve (9) for Sb. For the 
illustrative example shown in Fig. 2, where a = t'ext and b = K, (9) simplifies to 

Aiext'^^ext 

6K = ---=-, 

Ak V 

since S and T appearing in the respective A's cancel each other. 

To determine critical connections in a more complex model, we choose 6 = fc, i.e. the vector of 
indegrees, and solve (9) with a decomposition into eigenmodes. We can write the right-hand side as 


- (AkSk). = - ^ Tm Jy -I- TiJfj) VjSK-tj . 

3 


( 10 ) 


The equation holds because are only affected by connections to population i, and therefore their 

derivatives dfj-i/dki, dcji/dki and hence Ak^iu vanish for I ^[{i— 1)N + l,iN]. We now make the ansatz 


SKij — ^ 
i 


r„, (5J + T(J® J))y ’ 


( 11 ) 


which decomposes the changes SK into eigenmodes of the effective connectivity. The and are the 
l-th right and left eigenvectors of M as defined in (8), fulfilling the bi-orthogonality = Sin and the 

completeness relation = 1. Inserting (10) with (11) into (9) yields 


AaSa = — eiu^v^'^u. 
i 


(12) 
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Thus we can solve for ei by multiplying from the left with 


v’^^AaSa 
^ ^ > 

=:a'^ 


^ ^nl =:C'^ 

a" 

£,n ■ 


Our goal is to hnd a set of connections which dominate the size of the basin of attraction of the LA fixed 
point. Inserting (12) into (6) leads to 


Si^* 


I 

Sr-v 




St 


— A; 


-u 


where A; are the eigenvalues of M, which are either real or complex conjugate pairs since M £ . 

To determine the influence of each eigenmode on the shift of the fixed point, we project the eigenvectors 
onto the fixed-point shift Su* by multiplying each side with dv'*Si'*^ and solve again for du* to 
obtain 


Si^* 


—a 


8v* 




1 - Ai 5u*^5u* 


du*, 


-■m 


(13) 


where we define the (possibly complex-valued) coefficients fji. We aim at a decomposition of 5u* into 
real components. If Im(A/) = 0, fji is real, so we can work directly with rji := fji. Complex eigenvalues 
Im(Ai) 7 ^ 0 and corresponding eigenvectors come in conjugate pairs, so in this case we combine the 
corresponding coefficients rji := fji + fj*, to have all contributions rji € M. and = 1 by construction 

(13). Each r]i quantifies how much of the total fixed-point shift can be attributed to the l-th eigenmode, 
which allows identification of the most effective eigendirection (see Results), where we apply the ansatz 
(11) to a multi-area, multi-layer model of the vision-related areas of macaque cortex. 

The spiking simulations of the network model were carried out on the JUQUEEN supercomputer 
(Jiilich Supercomputing Centre, 2015). All simulations were performed with NEST version 2.8.0 (Eppler 
et al., 2015) with optimizations for the use on the supercomputer which will be included in a future 
release. The simulations use a time step of 0.1ms and exact integration for the subthreshold dynamics 
of the leaky integrate-and-fire neuron model (reviewed in (Plesser & Diesmann, 2009)). 
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Multi-area model 

The multi-area model of the vision-related areas of macaque cortex uses the microcircuit model of (Potjans 
& Diesmann, 2014) as a prototype for all 32 areas in the FV91 parcellation (Felleman & Van Essen, 1991) 
and customizes it based on experimental findings on cortical structure. From anatomical studies, it is 
known that cortical areas in the macaque monkey are heterogeneous in their laminar structure and can 
be roughly categorized into 8 different architectural types based on cell densities and laminar thicknesses. 
This distinction was originally developed for prefrontal areas (Barbas & Rempel-Clower, 1997), and then 
extended to the entire cortex (Hilgetag et ah, 2016). The visual cortex, and thus the model, comprises 
areas of categories 2, 4, 5, 6, 7 and 8. Precise layer-specific neuron densities are available for a number 
of areas, while for other areas, the neuron density is estimated based on their architectural type (see 
(Schmidt et ah, 2016) for details). 

The inter-areal connectivity is based on binary data from the CoCoMac database (Stephan et ah, 
2001; Bakker et ah, 2012; Felleman & Van Essen, 1991; Rockland & Pandya, 1979; Barnes & Pandya, 
1992) indicating the existence of connections, and quantitative data from (Markov et ah, 2014a). The 
latter are retrograde tracing data where connection strengths are quantified by the fraction of labeled 
neurons {FLN) in each source area. The original analysis of the experimental data was performed in the 
M132 atlas (Markov et ah, 2014a). Both the FV91 and the M132 parcellations have been registered to 
E99 space (Van Essen, 2002), a standard macaque cortical surface included with the software tool Caret 
(Van Essen et ah, 2001). This enables mapping between the two parcellations. 

On the target side, we use the exact coordinates of the injections to identify the equivalent area 
in the EV91 parcellation. To map the data on the source side from the M132 atlas to the EV91 par¬ 
cellation, we count the number of overlapping triangles on the F99 surface between any given pair of 
regions and distribute data proportionally to the amount of overlap using the F99 region overlap tool 
on http://cocomac.g-node.org. In the model, this FLN is mapped to the indegree Kab the target area 
A receives from source area B divided by its total indegree, i.e., FLNab = Kab/'^biKab'- Here, 
Kab is defined as the total number of synapses between A and B divided by the total number of neu¬ 
rons in A. On the source side, laminar connection patterns are based on CoCoMac (Felleman & Van 
Essen, 1991; Barnes & Pandya, 1992; Suzuki & Amaral, 1994; Morel & Bullier, 1990; Perkel et ah, 1986; 
Seltzer & Pandya, 1994) and on fractions of labeled neurons in the supragranular layers (SLN) (Markov 
et al., 2014b). Gaps in these data are bridged exploiting a sigmoidal relation between SLN and the 
logarithmized ratio of overall cell densities of the two areas, similar to (Beul et ah, 2015). We map 
the SLN to the ratio between the number of synapses originating in layer 2/3 and the total number of 
synapses between the two areas, assuming that only excitatory populations send inter-area connections, 
i.e., SLNab = Y^i Yi j ^ab^Y where the indices i and j go over the different populations 

within area A and B, respectively. In the context of the model, we use the terms FLN and SLN to refer 
to the relevant relative indegrees given here. On the target side, the CoCoMac database provides data 
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from anterograde tracing studies (Felleman & Van Essen, 1991; Rockland & Pandya, 1979; Jones et ah, 
1978; Seltzer & Pandya, 1991). 

Missing inputs in the model, i.e., from subcortical and non-visual cortical areas, are replaced by 
Poissonian spike trains, whose rate r'ext is a free, global parameter. In the original model all populations 
of a particular area receive the same indegree of external inputs AText- The only exception to this rule 
is area TH where the absence of granular layer 4 is compensated by an increase of the external input 
to populations 2/3E and 5E by 20%. To elevate the firing rates in the excitatory populations in layers 
5 and 6, we increase the external drive onto these populations. The possibility of a higher drive onto 
these populations is left open by the sparseness of the corresponding experimental data. We enhance the 
external Poisson drive of the 5E population, parametrized by the K^E,ext incoming connections per target 
neuron (indegree), in all areas by increasing k = AlsE.ext/Alext- The simultaneous increase in the drive of 
6E needs to be stronger, since the firing rates in population 6E of the original model (Pig. 3C) are even 
lower than the rates of 5E (averaged across areas: 0.09 spikes/s for 5E compared to 2 • 10“® spikes/s for 
6E). We thus scale up ATsE.ext linearly with k such that n = 1.15 results in ATsE.ext/Alext = 1-5. 
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